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The coupling parameter expansion in thermodynamic perturbation theory of simple fluids is gen- 
eralized to include the derivatives of bridge function with respect to coupling parameter. We applied 
the theory to Square- Well (SW) and Lennard-Jones(LJ) fluids using Sarkisov Bridge function. In 
both cases, the theory reproduced the radial distribution functions obtained from integral equation 
theory (IET) and simulations with good accuracy. Also, the method worked in the liquid-vapor 
coexistence region where the IETs are known to fail. In the case of SW fluids, the use of Carnahan- 
Starling expression for Helmholtz free energy density of Hard-Sphere reference system has improved 
the liquid-vapor phase diagram(LVPD) over that obtained from IET with the same bridge function. 
We also obtained the surface tension for SW fluids of various ranges. Results of present theory and 
simulations are in good agreement. In the case of LJ fluids, the equation of state obtained from the 
present method matched with that obtained from IET with negligible deviation . We also obtained 
LVPD of LJ fluid from virial and energy routes and found that there is slight inconsistency between 
the two routes. The two applications above lead to following conclusions. In cases where reference 
system properties are known accurately, the present method gives results which are very much im- 
proved over those obtained from the IET with the same bridge function. In cases where reference 
system data is not available, the method serves as an alternative way of solving the Ornstein-Zernike 
equation with a given closure relation with the advantage that solution can be obtained throughout 
the phase diagram with a proper choice of the reference system. 
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I. INTRODUCTION 



Describing accurately the thermodynamics and structure of fluids with short ranged inter-particle potentials has 
been a long standing problem in liquid state theory. Various methods have been developed to attack the problem. 
They can be divided into two categories. One category is based on integral equation theory(IET) and the other is 
based on thermodynamic perturbation theory (TPT) Methods based on IET take into account correlations between 
particles and aim at accurate description of structural properties of liquids i.e., the radial distribution function (RDF), 
direct correlation function (DCF) etc. The IETs involve solving the Ornztein-Zernike equation(OZE) coupled with 
a closure relation. Thermodynamic properties like internal energy and pressure are obtained as integrals of RDF. 
But inside the two phase region, most of the closures do not have a solution^. Apart from this, the IET method 
suffers from lack of thermodynamic consistency between the energy, virial and compressibility routes. Various closure 
relations were proposed to minimize the inconsistency @. Also closures with adjustable parameters have been proposed 
to enforce thermodynamic consistency by adjusting them. The requirement of consistency makes the computations 
complicated. Despite the improvements, accurate description of thermodynamic and structural properties of simple 
liquids with narrow potentials is still lacking. 

On the other hand, in earlier formulations of TPT, the structure and thermodynamic properties of reference system 
are assumed to be known and the thermodynamic properties of actual system are obtained as a small perturbation 
over the reference system. The structure of the reference system and the actual system are assumed to be the 
same. Within this approach, it was not possible to go beyond second order term in the perturbation series as they 
require higher order correlation functions [|[. ZhouQ developed a method in which the Helmholtz free energy of the 
system is written as a series over coupling parameter £. This method requires derivatives of the radial distribution 
function (RDF) at £ = 0. Zhou calculated the derivatives using finite-difference method and could obtain terms in 
the perturbation series up to fifth orderp. Beyond fifth order, it was difficult to obtain reliable estimates of the 
derivatives because of numerical complications. 

Recently, we attempted to address the problem in a more fundamental way combining the ideas of IET and TPT. 
We assumed that the RDF and the DCF of a system can be written as Taylor series in coupling parameter £ at £ = 0. 
We showed that the terms in Taylor series of both RDF and DCF (i.e., their derivatives) can be obtained by solving a 
simple set of equations obtained by differentiating the OZE and the generalized closure. The formula for free energy 
given by Zhou naturally emerges from our formalism apart from the RDF and DCF of the actual system. Within this 
approach we were able to obtain DCF, RDF and LVPDs for Square Well (SW) fluids of various ranges up to seventh 
order. By comparing the LVPDs obtained from third, fifth and seventh order versions of our theory, we concluded 
that the perturbation series for Helmholtz free energy has practically converged by seventh order. However, the series 
for RDF and DCF didn't converge for low densities and low temperatures for narrow square well fluids. Also, even 
though the Helmholtz free energy series has converged, we found that there was significant deviation of the obtained 
LVPD using seventh order TPT from simulation results for narrow SW fluids. The deviation and slow convergence 
could be because of two reasons: Firstly, we neglected the derivatives of bridge function w.r.t. £ in our calculations 
and secondly the error in the bridge function we have chosen for calculation. 

Also in our previous work, we were able to obtain the RDF and DCF even inside the LVPD and close to the 
critical region without any problems of convergence of the iteration scheme even though our numerical scheme is 
very simple. This raises questions about the existence of the solution to OZE inside the spinodal region. The RDF 
and DCF in the two-phase region find applications in the density functional theory of fluids @- One example is the 
calculation of surface tension. This calculation requires the direct correlation function in the two phase region. As 
the IETs doesn't have solution in some part of the LVPD, earlier calculations of surface tension were done using 
an interpolation formula given by Ebner et.ai.Q to obtain c(r) inside the two-phase region. In present paper, we 
address these issues with application to SW fluids. We use the bridge function proposed by Sarkisov[9j with slight 



modification for SW wellsfll] and include its derivatives w.r.t. C in the calculation. We calculate the RDFs for SW 
fluids of range 1.3 and LVPDs for SW fluids of ranges 1.25 and 1.375 using seventh order version of our method. 
Hard-sphere fluid is taken as reference system and Carnahan-starling[l^ expression is used for free energy density of 
the hard sphere reference system. The RDF and DCF of the reference system are obtained by solving the OZE within 
the Sarkisov approximation modified for SW wells [11]. Interestingly, we could obtain the RDF and DCF of the SW 
fluid even inside the spinodal region where the Sarkisov's closure is supposed to have no solution. We calculate the 
surface tension for SW fluids of ranges 1.25, 1.375 and 1.75 using the expression obtained from the square gradient 
functional for Helmholtz free energy and compare with simulation results. 

As an application of our theory to non-hard sphere reference systems, we apply it to Lennard- Jones (LJ) fluid. We 
use the Sarkisov bridge function for both reference system and perturbation part. The RDFs, isotherms and LVPD 
for LJ fluid are obtained and compared with those obtained from IET and simulations wherever available. The paper 
is organized as follows. In Section II we discuss the method briefly. In Section III we apply the theory to SW fluids 
and in Section IV to LJ fluids and the results are analyzed. The paper is concluded in Section V. 
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II. THEORY 



Consider a fictitious system at temperature T and density p with interaction potential it(£, r) given by 

u(C,r) = u ref (r) +(u pert (r) (1) 

£ is a coupling parameter. When £ = 0, the fictitious system becomes reference system with interaction potential 
u re f(r). A non-zero C will add a perturbation to u re f(r) as shown in above equation. £ = 1 corresponds to the 
potential of the system under consideration. We postulated that the RDF and DCF of the system with potential 
r) can be written as a McLaurin series in that is, 

dc \ C 2 fd 2 c 



^r) = c o( r) + C[ M ) ( __ o + ^ [w ) c __ Q + (2) 

KC,r)=*(r) + C('^ + + (3) 



Hereafter we shall denote n th partial derivative of any function r) w.r.t Q as X r). With simple arguments, 
it can be seen that the Taylor series expansion for c(£, r) and g{C,: r ) w.r.t. £ exists. This requires to prove that both 
c(C, r) and </(£, r) are analytic in £. Using diagrammatic techniques, DCF in general can be written as follows [l[ 

c(l,2) = {sum of all topologically distinct diagrams that consist of two white 1-circles labelled 1 and 2, (4) 
black -circles and fu -bonds and are free of connecting circles} 

Notation and conventions are followed as given in ref.([J|). Where pS^ is the single particle density which is equal to 
density p in a homogeneous system. /m(U 2) = e - ' 3 "' 1,2 ' — 1 is the Mayer's function. u(l, 2) is the potential between 
the particle 1 and particle 2 defined by Eq.([T]) (£ is not shown for notational convenience). /3 = l/fc^T. Clearly, 
potential u(l,2) enters the definition of c(l,2) through the Mayer's function in each term. Separating the reference 
and perturbation part of the potential in the exponential term of the Mayer's function, 

/m(1, 2) = e (-£(«^(i,2))) e (-/3(CWi,2))) _ ! (5) 

From the above Eq.([5]), it can be seen that the Mayer function is analytic in £ unless u pert (l, 2) has some singularity, as 
the exponential function is an analytic function. Thus the Mayer's function, which is the first term in the diagrammatic 
expansion of c(l,2) can be written as a Taylor series in £ which is nothing but the series obtained by expanding the 
exponential containing £. 

Since the Mayer's function is the basic component in all the diagrams of c(l,2), all other diagrams also can be 
expanded in powers of £ and thus Taylor series expansion exists for each diagram and hence c(l,2) can be expanded 
as a Taylor series in (. In the case of homogeneous fluid, c(l,2) depends only on the distance between particles and 
is written as c(£, r) where r is the distance between the particles. 

Also, the diagrammatic expansion for RDF is 

g(l,2) = {sum of all topologically distinct diagrams that consist of two white 1-circles labelled 1 and 2, (6) 
black -circles and fu-bonds and are free of articulation circles] 

From similar arguments given in the case of DCF, it can be seen that g(Cj r ) is analytic in Q and a Taylor series 
expansion exists. Now, with the understanding that a Taylor series expansion exists for both c((,r) and g((,r), we 
proceed to obtain relation between their derivatives as follows: Using the concept of potential of mean force, the 
radial distribution function can be written as 

9 (C,r) = exp(</>(C,r)) (7) 
0(0 r) = -0 (u ref {r) + (u pert (r)) + y(£ r) + B(£ r) 

where r) is the indirect correlation function defined as h((,r) — c(£, r) and h((,r) — g(C> r ) — 1 is the total 
correlation function of the fictitious fluid. The bridge function B(£,r) is a sum of an infinite series of the "bridge 
diagrams" [3[. Since g(C, r )> and hence h(£,r), as well as c(£,r) are expanded in a series in £, the correlation function 
y(C, r) is also a series in £. The n th order coefficient in its series is given by y^(Ci r ) — (C r) — a 1 -' 1 ' (£, r). Generally, 
r) is chosen as a function of r). Several approximations to J5(^,r) in terms of y(Ci r ) an d certain empirical 
parameters are available 
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The general expression for n th order derivative of g{(,, r) from Eq.® is 

(n-l) 

9 {n HCr) = £ \ C ^ 1] \ (n - m) (C,r) . 9 (m) (0), n > 1 (8) 

m=0 

where [Cm ^] is the binomial coefficient. The derivatives (f>^ n \Ct r ) are given by 

^")(C,r) = - y 8u P ert(r)<5„,i+j/ (n) (C,r) + J B(")(C,r), » > 1 (9) 

where <5 nj i is the Kronecker delta. The derivatives g^ n \C> r ) can be computed using Eq.® in a recursive manner, 

using 0( n ) (£, r) either from initial guess or previous iteration. 

Since 

c(Cr)=g((,r)-y((,r)-l, (10) 

its n th order derivative is 

c^(Cr)=g^(C,r)-y^(C,r), n>l. (11) 

To get another set of relations between cW(£, r) and j/"* 1 ^, r), we consider the Ornstein Zernike Equation (OZE) 
in Fourier space: 

hK - k > = T^h) (12 » 

where k) and c(£, k) are the Fourier transforms of h(Q, r) and c(C, r), respectively. 
Also, structure factor s((, k) is defined as 

Therefore, using above equation in OZE i.e., Eq.(fT2j) we get 

s((,k) = l + ph((,k) (14) 

Now, differentiating Eq. (fT3]) n-times, we get 

(n-i) 

S (")(C,fc) = [s (O) (C,£0] p ^- m >(C,fc) S (m) (C,fc), n> 1 (15) 

m— 

Using Eq.lpHp the derivatives of y(C, A;) = k) — c(C, fc) are expressed as 

2/ W(C,fe) = p- 1 S W(C,fe)-cW(C I fe), n>l (16) 

Thus in order to obtain up to N th order term in Taylor series expansion of both c(l, r) and r), TV coupled linear 
equations in real-space and AT coupled linear equations in Fourier space have to be solved simultaneously with £ = 0. 
For example, to obtain Taylor series expansion up to second term, the set of four equations given by 

c^{0,r) = (-pu per t(r)+y il) (0,r)+B^(0,r))g(0,r)-y^(0,r) (17) 

c^(0,r) = (-/5v*W + !/ (1 \0^)+5 ( ^0, r ))VO, r ) + B (2) (0 I r)j(0,r) + ! / l 2 \0,r)( S (0,0 - 1) (18) 

y (1) (0, k) = c (1) (0, fc)(s 2 (0, k) - 1) (19) 

y (2) (0,fc) = c( 2 >(0,fc)(s 2 (0,fc) - 1) + 2p(c«(0, k)) 2 s 3 (0,k) (20) 

have to be solved (where g(0, r) is y(0, r) + c(0, r) + 1 ). 

The CPE for Helmholtz free energy density f{T, p) of a homogeneous fluid at a given temperature T is given by [l[ 

p 2 n 



/(T,p) = / re/ (p) + ^ / dC / rfr Upert (r).g(C,r) (21) 
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where f re f (p) is the free energy density of the reference system. For notational simplicity we do not show temperature 
dependence of f(T,p) explicitly hereafter. Substituting Eq.Q in Eq. (|2"T) and integrating over £, we get 

f(p) = f ref (p) + dru pert (r) L(r) + ± 5 (1) (0,r) + I </ 2 )(G,r) + \ (22) 

Here we have used the shortened notation for the derivatives (d n g(r)/d( n ) ( * =0 which is readily obtained as y^(0, r) + 

c ( n )(0,r). Thus the method provides the DCF, RDF as well as the free energy density. The Taylor series of RDF 
truncated up to second order gives third order CPE for f(p). 

A. Numerical Procedure 

The numerical method is simple and is as explained in our previous paper. However, for the sake of completeness 
we explain it here. Firstly, it has to be noted that the present method offers the flexibility to use any known quantity 
of the reference system like c'°'(0,r) = Co(r), g^(0,r) = go( r ) or the reference system free energy density f re f(p). 
For example, in the case of Hard-Sphere reference system, the Carnahan-Starling expression for f re f(p) c&n be used. 
Similarly, <?o( r ) and Co(r) obtained from simulations may be used. Alternatively, they can be generated by solving the 
OZE using some bridge function using the same procedure as outlined below. In such a case, the method amounts 
to solving the OZE in an alternative way. To compute the derivatives, say, up to n th order, Eq. lTTTj) in real space 
and Eq. (|16[) in Fourier space are solved employing an iterative procedure. First of all, we choose guess solutions 
for y( m \0,r) (in practice a null solution suffices), for all m in the range 1 < m < n, and compute g( m ^(0,r), using 
Eq.©. Then, C ( m '(0,r) are obtained using Eq. (fTTj) . Next, their Fourier transforms c^^Ojfc) are computed using 
an FFT algorithm. Mesh widths in the range 0.01 to 0.001 are found to adequate. These are employed in Eq. (TT5")) 
to obtain s^ m ^(0, k). These functions, when used in Eq. tfTB")) . provide y( m \0, k). Inverse Fourier transformation of it 
gives ?/ m )(0,r). This completes the first iteration. However, before starting the next iteration, we employ a linear 
mixing of the initial and new solutions: ay( m \Q,r) + (1 — a)y^ m ' (0, r) — > i/ OT )(0,r). They are then used in Eq. fTTj) 
for the second iteration. The procedure is repeated until the root mean square differences between successive iterates 
of y (™)(0,r), for all m, are less than a prescribed tolerance. We find that a — 0.5 and tolerance of 10 8 are adequate 
to get accurate solutions. 

III. APPLICATION TO SQUARE WELL FLUIDS 

Seventh order version of above mentioned theory is applied to SW fluids. We use the closure proposed by Sarkisov 
with slight modification by MendoubflTl]. 

5(0 r) = (1 + 2y*(C,r)) 1 / 2 - 1 -y*((,r) (23) 

where 

y*(C.r) = y(C,r) + pf M (a+)/2, r < a (24) 
= y(t,r)+pf M {r)/2, r>a 

and a is the Hard sphere diameter. We obtained the RDF and DCF of the reference system using the same bridge 
function by solving the OZE using a similar procedure as explained above. 

Reduced units (e/ks = u = 1, where e is the well depth and a is hard sphere diameter) are used throughout the 
paper. In Fig.© we compare g(r) obtained using 7 th order version of our TPT and that obtained through lETpjJ 
for SW fluid of range 1.3 for densities 0.2 and 0.8 and temperature T — 1.0. Clearly, there is negligible difference 
between results obtained using present method and those obtained from IET. We observed that except for very low 
temperatures and low densities, results of fifth order and seventh order TPTs have negligible deviation showing that 
the Taylor series has converged. Also convergence of the iteration scheme is good in the whole phase diagram. g(r) 
obtained by our method for SW fluid of range 1.25 in the spinodal region at p = 0.4 and T — 0.65 is shown in Fig.([T|). 
Above observations show that the present method can be viewed as an alternative way of solving the OZE. Advantages 
of solving the OZE using the perturbation method are, firstly solutions to OZE can be obtained inside the spinodal 
region also where earlier methods fail and secondly a simpler numerical scheme is sufficient throughout the phase 
diagram. Apart from these, any known information about the reference system can be used to improve the accuracy 
of the calculations. For example, in the present case we used the Carnahan-Starling expression for Helmholtz free 
energy density of the reference system. In Fig.© and Fig.© we give plots of LVPD for SW fluids of ranges 1.375 
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and 1.25 respectively using seventh order version of TPT. Our results are compared with simulations, those obtained 
from IET[11] and our previous results using Malijevsky-Labik brdige function neglecting the derivatives of bridge 
function. Results of the present calculations are very much improved in the liquid part of the LVPD over those of 
Mendoub. The main reason for this improvement is the use of Carnahan-Starling expression for free energy. Also the 
present results significantly improved our previous results [12] where we neglected the derivatives of bridge function. 
Possibility of obtaining c(r) in two phase region opens up applications in density functional theory(DFT) of liquids. 
One simple application is studying the surface properties of liquids like surface tension etc. We use the c(r) obtained 
from seventh order version of TPT to obtain surface tension for SW fluids. Formula for surface tension is obtained 
by Yang et. al. [l3T ] from square-gradient functional for Helmholtz free energy of inhomogeneous fluids. A brief 
derivation is as follows: The square-gradient approximation for Helmholtz free energy functional of an inhomogeneous 
liquid occupying volume f2 at temperature T is given by 



F[p(r)} = J dr {/(p(r)) + / 9 (p(r))|Vp(r)| 2 } (25) 

where f(p) is Helmholtz free energy density of homogeneous liquid. The second term is the effect of inhomogeneity. 
p(r) is the number density in an infinitesimal volume around r.F[p(f)], f(p(r)) and fg(p(r)) are all functions of T 
even though the dependence is not shown explicitly. 

The coefficient f g of gradient term also referred as influence parameter is 

fM?)} = ^ J d^r"c[p(f),r'} (26) 

We assume that liquid-vapor interface is flat and that z-axis is normal to the interface pointing out into the vapor 
from the liquid. In such a case, Eq. (|27|) becomes 



F[p(z)]=Aj dz {/(p(z)) + / g (p(z))|^p } 

where A is the surface area. Grand free energy of the system is given as 

r[p(z)} = F[p(z)} - fiN 

where p is the chemical potential of the system and N is the total number of particles 
Minimizing T w.r.t. p we get 



d_ 

dz 



fM 



dz 



where j(p) is the grand free energy density. 

Integrating above Eq. (j2"9")l with boundary conditions 



p(z — > oo) = p g ; p(z — > — oo) = pi and 



dp(z — > ±oo) 
dz 



gives 



dp 

dz 



l(p) - l(Pl) 



fM 



1/2 



(27) 



(28) 



(29) 



(30) 



(31) 



where pi and p g are coexisting liquid and vapor densities for temperature T under consideration. Surface tension S 
can be calculated using the formula 



Using Eq. (|3"Tj) in above equation gives 



5 = 2 £>(£)"* 



S = 2 [f g (p)(l( P )-7(pi))} 



(32) 



(33) 
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In the above equations j(pi) may be replaced by j(p g ) also as both have same value. 

Above explained formalism is applied to SW fluids of ranges 1.375, 1.5 and 1.75 with c(r) obtained from seventh 
order TPT as explained above. Fig. (J4j) shows f g (p) as a function of p for different temperatures for SW fluid of range 
1.25. From Fig.(Q} it can be seen that f g (p) becomes negative at high densities which is unphysical. This might be an 
artifact of the approximate bridge function used. However, within the binodal, where f g is required for the calculation 
of surface tension(S), it is positive. Surface tensions obtained from our calculation arc plotted as a function of T in 
Fig.©. Our results compare well with simulations except close to the critical temperatures. 

IV. APPLICATION TO LENNARD JONES FLUID 

The theoretical formalism explained in section II is general and can be applied to non-hard sphere reference systems 
also. As an example, we apply the theory to Lennard Jones(LJ) fluid. The Lennard Jones potential is split into 
reference and perturbation parts according to the Weeks Chandler and Anderson (WCA) [lij method. For this case, 
we use for -B(C, f ), the expression proposed by Sarkisov Q which is same as Eq.(|23j) with 

V*(C, r ) = y((, r ) ~ pf3u{r m ), r <r m (34) 
= y(C,r) - p(3u{r), r>r m 

where r m is the minimum energy point of the Lennard- Jones potential. The reference system c(r) and g(r) are 
obtained by solving the OZE with the same bridge function. Thus, in effect, we solve the OZE for LJ fluid with the 
above bridge function using the perturbation method. In Fig.©, we compare the RDF of LJ fluid obtained from 
seventh order TPT with simulation results. There is excellent agreement between seventh order TPT and simulation 
results for the cases shown except for a slight deviation for the case with lowest temperature. In Fig.©, Equation of 
State(EOS) of LJ fluid for various temperatures obtained using seventh order TPT is compared with those obtained 
by solving OZE by Sarkisov[9j. Pressure(P) is obtained using the virial formula given by 

P = pk B T - l -p 2 J™ ^lg( r )4nr 2 dr (35) 

Values obtained by our method matched with negligible deviation from those given by Sarkisov. Also, using our 
method we could obtain the pressure(P) for all density points in the phase diagram as g(r) at any point in the phase 
diagram could be calculated. This is an advantage of the present method over directly solving the OZE as solution 
cannot be found inside the two phase region. Again, even close to the critical region, the same numerical procedure 
was sufficient whereas IETs face numerical convergence problems in this region. In Fig.®, LVPD of LJ fluid is shown 
and compared with simulation results. We plotted the LVPDs obtained in two ways. One is obtained by Maxwell 
construction of pressure isotherm obtained using Eq. ([35|) which is the so called virial route. This is shown in solid 
line in Fig.®. LVPD shown in dashed line in Fig.® is obtained using the energy route. It is obtained as follows: 
Pressure isotherm of the reference fluid is obtained using Eq.p5p with u re f(r). From this, Hclmholtz free energy 
density of the reference fluid is obtained from the formula below 

fref(p) = P k B T l\ Pref i pl) - 1) + P k B Tln(p) (36) 
Jo P 

Once the reference free energy is obtained, method described in section II can be applied to get the free energy density 
of the required system. The pressure isotherm is obtained by differentiating the Helmholtz free energy w.r.t. volume 
of the system. Maxwell construction is done to get the coexistence points. From the figure, it can be seen that the 
LVPDs obtained from both the routes differ slightly along the liquid part of coexistence curve. Also, there is some 
deviation of both LVPDs from the simulation results. This is due to the bridge function used. Even though Sarkisov 
bridge function is supposed to be quiet accurate, there is still some inconsistency between various thermodynamic 
routes. Imposition of the thermodynamic consistency between various routes as a constraint would solve the problem 
and may improve the accuracy of the results. This can be done in a straight-forward way within the theory presented 
in the paper. 

V. CONCLUSIONS 

In this paper we investigated in detail the coupling parameter expansion of simple fluids. We have generalized our 
previous version by including the derivatives of the bridge function in the method. The theory has been applied to 
study square well fluids and Lennard Jones fluids using Sarkisov bridge function. 
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In the case of SW fluids, we showed that the g(r) obtained from seventh order TPT matches well with that obtained 
from IET. We noticed that solutions to OZE can be obtained inside the spinodal also where ordinary IETs fail. The 
LVPDs obtained using the present method improved over those obtained from IET. This is because of usage of 
accurate parametrization of reference system free energy density which is not possible in IET. Usage of a improved 
bridge function and including the derivatives of the bridge function also in the coupling parameter expansion improved 
the LVPDs of SW fluids over our earlier calculations. Possibility of obtaining solutions of OZE in two phase region 
helped to obtain interfacial properties. We obtained the surface tension of SW fluids of ranges 1.375, 1.5 and 1.75. 
Agreement with simulation results is found to be quiet good except close to the critical region. As an application to 
non hard-sphere fluids, we applied the method to Lennard- Jones fluids. In this case also, we found that the seventh 
order perturbation theory results match with those obtained using IET with the same bridge function. Additionally, 
we were able to obtain solutions in the spinodal region also where the IET doesn't have solutions. LVPD is obtained 
through energy and virial routes and compared with simulations. 

Thus present method can be viewed as an alternative way of solving the Ornstcin Zernike equation with a given 
closure. The advantages that we noticed in solving the OZE in this way are, firstly, the solution can be obtained in the 
coexistence region also where earlier methods of solution fail. Secondly, a simple numerical scheme is sufficient and 
works even in the critical region. Apart from this, the perturbation method allows to choose a reference system whose 
properties may be known more accurately. Also, the reference system can always be chosen in such a way that its 
interaction potential u re f(r) is only repulsive and thus doesn't have a liquid vapor coexistence. For such a reference 
system, any closure would give solution in the entire phase diagram. Since the perturbation method described in the 
paper confines to obtaining derivatives strictly at £ = 0, and the n th order derivative depends only on derivatives 
upto (n — l)*' 1 order, all the derivatives exist once the reference system has a solution throughout the phase diagram. 
More clearly, the first derivatives of RDF and DCF of the required system depend only on the reference system 
RDF and DCF. Hence they can be obtained if reference system RDF and DCF are known. Once they are obtained, 
second derivatives of RDF and DCF can be calculated as they depend only on reference system RDF, DCF and first 
derivatives of RDF, DCF and so on. Thus the existence of solution of OZE for the actual system depends on the 
existence of solution of OZE for the reference system which will in general exist if the reference system is chosen to 
have only repulsive interaction. The convergence of Taylor series expansion is seen to be quiet satisfactory except for 
very low temperatures and low densities for narrow wells. Thus the accuracy of the present method is limited only 
by the bridge function used. Implementation of consistency condition in bridge function may improve the accuracy 
of the present method which is a subject of our future work. 
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FIG. 1. (Color OnlinejTop: g(r) for SW fluid of range 1.3 for densities 0.2,0.8 and T = 1.0. (Circles: simulation results [HI ); 
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T = 0.65 and density 0.4. 
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FIG. 2. (Color Online) LVPD of SW fluid of range 1.375. ( dotted line: 7th order TPT with Malijevsky Labik bridge 
function [12]); (solid line: 7th order TPT with Sarkisov B(r) (present work)); (Short dashes IET results with Sarkisov B(r) 
obtained from Mendoub et. al.[ll]]);( Squares and stars are simulation results [la. fl7j): 
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FIG. 3. (Color Online)LVPD of SW fluid of range 1.25. ( dotted line: 7th order TPT with Malijevsky Labik bridge function[H]); 
(solid line: 7th order TPT with Sarkisov B(r) (present work).); (Short dashes IET results with Sarkisov B(r) obtained from 
Mendoub et. al.);( Squares and stars are simulation results [lq. Il7fl): 



11 




1.2 



\ 




0.6 0.8 1 1.2 1.4 1.6 1.8 



Temperature 

FIG. 5. (Color Online) Surface tension for SW fluids in reduced units. From left to right are for ranges 1.375, 1.5 and 1.75. 
stars, pluses, crosses are simulation results [l^.[l9|. 



12 




r r 



FIG. 6. (Color Online) LJ rdfs for different temperatures and densities obtained using 7 th order TPT. (pluses: Simulation 
results[lj), (Top left: T = 1.552, p = 0.45), (Bottom left: T = 2.934, p = 0.45), (Top right: T = 0.658, p = 0.85), (Bottom right: 
T = 2.888, p = 0.85). T and p are in reduced units 
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FIG. 8. (Color Online) LVPD of LJ fluid. Dashed line is obtained from EOS calculated using energy route. Solid line is from 
EOS obtained using virial route. In both cases Sarkisov bridge function and 7 th order TPT are used. Crosses are simulation 
results [H]]. 



